Conducting Meta-Analytic Structural Equation Modeling with R

My background

• Academic position:
  • Associate Professor at Department of Psychology, the National University of Singapore (NUS),
  • Associate Professor (courtesy appointment) at Department of Management & Organisation, NUS.
• Research interests:
  • Quantitative methods,
  • Integrating meta-analysis into the structural equation modeling (SEM) framework.
• Associate editors:
  • Research Synthesis Methods
  • Neuropsychology Review
  • Frontiers in Psychology (Quantitative Psychology and Measurement)
• Editorial boards:
  • Psychological Methods
  • Psychological Bulletin
  • Health Psychology Review
  • Journal of Management
• Special issues:
  • Co-edited a special issue with Adam Hafdahl on Meta-Analytic Structural Equation Modeling appeared in Research Synthesis Methods.
  • Co-editing a special issue with Suzanne Jak on Big Data in Psychology that will appear in Zeitschrift fuer Psychologie

Proposed schedule

• Introduction of MASEM
• Fixed-effects model
• Random-effects model
• Conclusions
• Q&A

R packages

• R project
• RStudio: A powerful IDE for R. It is useful, but not compulsory, to have it for this workshop.
• The metaSEM package
  • Website
  • Examples
  • Paper
• All workshop materials are available at https://github.com/mikewlcheung/masemWorkshop.
• The analyses were based on R 3.5.1, OpenMx 2.11.5, and metaSEM 1.2.0.
• You may need to install the following packages in R before you can run the analyses.

## Required packages for this workshop
lib2install <- c("metaSEM", "semPlot", "readxl")

## Install them automatically if they are not available on your computer
for (i in lib2install) {
  if (!(i %in% rownames(installed.packages()))) install.packages(i)
}

Some resources

Suggested readings

• Fixed-effects MASEM
• Fixed- and random-effects MASEM
• Alternative models for random-effects MASEM
• Applications and issues
• Some reflections on combining meta-analysis and SEM

Getting help in using MASEM

• metaSEM in OpenMx forum

What is structural equation modeling (SEM)?

• SEM represents a family of related techniques.
• Many popular multivariate techniques are special cases of SEM.
• It is also known as covariance structural analysis, covariance structure model, analysis of covariance structures, analysis of correlation structure, LISREL model (in the old days), etc.
• Both Karl Joreskog and Peter Bentler received the American Psychological Association Award for Distinguished Scientific Applications of Psychology at 2007:

“for their development of models, statistical procedures, and a computer algorithm (LISREL by Karl G. Joreskog; EQS by Peter M. Bentler) for Structural Equation Modeling (SEM) that changed the way in which inferences are made from observational data by permitting hypotheses derived from theory to be tested.”

Problems in primary research using SEM

• Low statistical power in psychological studies;
• Contradictory findings with significant tests;
• Confirmation bias: reluctant to consider alternative models in SEM;
• Different models proposed by various researchers.
• Conducting more empirical studies may not necessarily decrease the uncertainty of a particular topic if the findings are inconsistent (National Research Council, 1992).¹
• MASEM is used to address these issues.

Terms used in the literature

• Several names have been used interchangeably in the literature,
  • meta-analytic path analysis;
  • meta-analysis of factor analysis;
  • meta-analytical structural equations analysis;
  • path analysis of meta-analytically derived correlation matrices;
  • SEM of a meta-analytic correlation matrix;
  • path analysis based on meta-analytic findings; and
  • model-based meta-analysis.
• We use the generic term meta-analytic structural equation modeling (MASEM) to describe this class of techniques.

Brown and Stayman (1992)

• Topic: Antecedents and consequences of attitudes toward the advertisement (Ad)²
  • No. of variables: 5
  • No. of studies: 47
  • Pooled sample size across studies: +4,600

Premack and Hunter (1988)

• Topic: Individual unionization decisions³
  • No. of variables: 6
  • No. of studies: 14
  • Pooled sample size across studies: +2,800

Norton et al. (2013)

• There are ten different factor structures on the Hospital Anxiety and Depression Scale supported by some empirical data.⁴
• The authors identified 28 independent samples from 21 studies (N=21,820).
• They found that the bifactor structure consisting of a general distress factor and anxiety and depression group factors fitted the data best.

Murayama and Elliot (2012)

• Based on their meta-analysis, these authors found that the association between competition and performance was close to zero.⁵
• They proposed two mediators (performance-approach goals and performance-avoidance goals) to explain this apparent zero correlation.
• The total effect from competition to performance is close to zero because the indirect effects were positive and negative (474 studies with 139,464 participants).

Approaches to MASEM

• Nearly all methods use a two-stage approach to conducting MASEM:
  • Stage 1 analysis: Combine the correlation matrices into a pooled correlation matrix;
  • Stage 2 analysis: Fit structural equation models on the pooled correlation matrix.
• Common approaches:
  • Univariate approach (Hunter & Schmidt, 1990; Viswesvaran & Ones 1995);⁶
  • Generalized least squares (GLS; Becker, 1992);⁷
  • Two-stage structural equation modeling (TSSEM; Cheung, 2014; Cheung & Chan, 2005).⁸, ⁹
• Cheung (2018) provides personal (probably biased) reflections on combining meta-analysis with SEM.¹⁰

Univariate approach (Viswesvaran & Ones, 1995)

• Stage 1 analysis
  • It synthesizes the correlation coefficients in a correlation matrix as if the correlation coefficients were independent.
  • Incomplete correlation coefficients are handled by either listwise deletion or pairwise deletion (more popular).
  • Problem: when pairwise deletion is used, missing correlations are assumed missing completely at random (MCAR).
• Stage 2 analysis
  • The pooled correlation matrix is used as if it was a covariance matrix.
  • Researchers usually use the harmonic mean of the sample sizes as the sample size.
• Problems:
  • The correlation matrix is analyzed as if it was a covariance matrix. The test statistics and the SEs may be incorrect (Cudeck, 1989).¹¹
  • A single sample size is used to represent the precision of the average correlation matrix.
• The univariate approaches are based on the assumption that MASEM=MA+SEM.

One (sample) size does not fit all

• There are some problems with fitting SEMs with a single sample size in the univariate approach.
• Let us consider the average correlation matrix on the right as an example.
• The sample sizes for the correlations \(r_{21}\), \(r_{31}\), and \(r_{32}\) are 200, 500, and 1,000, respectively.
• The arithmetic mean, the harmonic mean, and the median are 567, 375, and 500, respectively.
• If the harmonic mean (375) is used as the sample size in SEM, some estimated SEs are over-estimated while others are under-estimated.
• The sample size affects the chi-square test, some goodness-of-fit indices, and SEs.

Problems of treating a correlation matrix as a covariance matrix

• Analysis of the covariance matrix:
  • No. of pieces of information: \(4 \times (4+1)/2=10\); no. of the unknown: 9, and df: 1 (10-9)
  • \(\hat{\Sigma} = \hat{\Lambda} \hat{\Phi} \hat{\Lambda}^T + \hat{\Psi}^T\)
  • The diagonals of \(\hat{\Sigma}\) are close but not necessarily equal to the variances of the variables.
• Analysis of the correlation matrix:
  • The diagonals are always 1. They do not carry any information.
  • No. of pieces of information: \(4 \times (4-1)/2=6\); no. of the unknown: 5, and df: 1=6-5.
  • The error variances are not parameters.
  • The diagonals of \(\hat{\Sigma}\) are precisely 1.
• The chi-square test or SEs can be incorrect if we treat a correlation matrix as if it was a covariance matrix in SEM.
• Cudeck (1989) warned the problems of treating correlation matrices as covariance matrices in SEM almost 30 years ago. However, many researchers still ignore the warnings in MASEM today.

GLS approach (Becker, 1992)

• Stage 1 analysis
  • GLS is used to synthesize correlation matrices by taking the dependence of the sampling error into account;
• Model:
  • \(r = X \rho_\mathrm{F} + e\), where \(r = \begin{bmatrix} r_1 \\ r_2 \\ \vdots \\ r_k \end{bmatrix}\), \(X = \begin{bmatrix} X_1 \\ X_2 \\ \vdots \\ X_k \end{bmatrix}\) and \(e = \begin{bmatrix} e_1 \\ e_2 \\ \vdots \\ e_k \end{bmatrix}\) with \(V=\begin{bmatrix} V_{1} \\ 0 & V_{2} \\ \vdots & \ddots & \ddots \\ 0 & \cdots & 0 & V_{k} \end{bmatrix}\).
  • The estimated correlation vector with its sampling covariance matrix are:
  • \(\begin{split} \hat{\rho}_\mathrm{F} & = (X^TV^{-1}X)^{-1} X^TV^{-1}r, \\ \hat{V}_\mathrm{F} & = (X^TV^{-1}X)^{-1}. \end{split}\)
• If a random-effects model is used, a different \(V^*\) which is usually larger than \(V\) is required.

GLS approach (2)

• Stage 2 analysis
  • The parameter estimates are obtained by matrix calculations whereas the SEs are calculated by the multivariate delta method.
  • Assuming that the first variable is the dependent variable and that the remaining variables are the predictors, the standardized regression coefficients in a regression analysis can be estimated by
  • \(\begin{split} R & = \begin{bmatrix} 1 & \\ R_{01} & R_{11} \end{bmatrix} \\ \hat{\beta} & = R_{11}^{-1} R_{01} \end{split}\).
• Practical issue:
  • It is not straightforward to generalize it to over-identified models or models with latent variables.

TSSEM approach

• The TSSEM approach is very similar to the GLS approach, but it is based on SEM.
• Fixed-effects TSSEM (Cheung & Chan, 2005)
  • Researchers are only interested in the studies included in the meta-analysis;
  • A common correlation/covariance matrix is assumed.
• Random-effects TSSEM (Cheung, 2014)
  • Studies are samples of a larger population;
  • Researchers may want to generalize the findings beyond the studies included;
  • Studies may have their own correlation/covariance matrices.

Fixed-effects TSSEM

• The distribution theory of SEM is based on covariance matrices, whereas correlation matrices are usually used in MASEM.
• We create additional variables (D) so that SEM can correctly analyze correlation matrices.
  • \(\Sigma (\theta) = DP(\theta)D\),
  • where \(\Sigma(\theta)\) is the structural model on the covariance matrix, \(D\) is a diagonal matrix, and \(P(\theta)\) is the structural model on the correlation matrix with the constraints that \(Diag(P(\theta)) = 1\) with \(1\) is a vector of ones (Joreskog & Sorbom, 1996).¹²
• For example, \(\Sigma (\theta) =\begin{bmatrix} 4.0 \\ 1.8 & 9.0 \\ 3.2 & 6.0 & 16.0\end{bmatrix}=\begin{bmatrix} 2 \\ 0 & 3 \\ 0 & 0 & 4\end{bmatrix} \begin{bmatrix} 1 \\ 0.3 & 1 \\ 0.4 & 0.5 & 1\end{bmatrix}\begin{bmatrix} 2 \\ 0 & 3 \\ 0 & 0 & 4\end{bmatrix}\).
• \(D\) is required to make sure that the distribution theory of covariance applies to correlation matrices.

Fixed-effects TSSEM: Stage 1 analysis (1)

• Analysis of correlation matrices with a multiple-group SEM
  • Common correlation matrix: \(P_\mathrm{F}=P_1=P_2=\ldots=P_k\), where \(D_i\) may vary across studies.
  • Test of homogeneity of correlation matrices: \(H_0: P_1=P_2=\ldots=P_k\) compared to a model that \(P_i\) may vary across studies.
  • Incomplete or missing correlations are efficiently handled by maximum likelihood estimation (Muthenn, Kaplan, & Hollis, 1987).¹³
  • The assumption of the homogeneity of correlation matrices may also be tested by the goodness-of-fit indices, such as the RMSEA and SRMR.

Fixed-effects TSSEM: Stage 1 analysis (2)

• Analysis of covariance matrices (Cheung & Chan, 2009)¹⁴
  • Common covariance matrix: \(\Sigma_\mathrm{F}=\Sigma_1=\Sigma_2=\ldots=\Sigma_k\)
  • Test of homogeneity of covariance matrices: \(H_0: \Sigma_1=P_2=\ldots=\Sigma_k\) compared to a model that \(\Sigma_i\) may vary across studies.

Fixed-effects TSSEM: Stage 2 analysis (1)

• After the first-stage analysis, the estimated common correlation matrix \(R_\mathrm{F}\) and its asymptotic sampling covariance matrix \(V_\mathrm{F}\) are available.
• When handling square matrices, it is easier to convert them into vectors, e.g.,
  • \(X = \begin{bmatrix}1 & 2 & 3 & 4\\ 2 & 5 & 6 & 7\\ 3 & 6 & 8 & 9\\ 4 & 7 & 9 & 10 \end{bmatrix}\),
  • \(vech(X)=\begin{bmatrix}1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10\end{bmatrix}^T\) for covariance matrix, and
  • \(vechs(X)=\begin{bmatrix}2 & 3 & 4 & 6 & 7 & 9\end{bmatrix}^T\) for correlation matrix.

Fixed-effects TSSEM: Stage 2 analysis (2)

• Suppose the proposed correlation structure is \(\rho_\mathrm{F}(\theta)=vechs(P_\mathrm{F}(\theta))\), the discrepancy function based on the weighted least squares (WLS) is (Bentler & Savalei, 2010):¹⁵
  • \(F_\mathrm{WLS}(\theta) = (r_\mathrm{F} - \rho_\mathrm{F}(\theta))^T V_\mathrm{F}^{-1} (r_\mathrm{F} - \rho_\mathrm{F}(\theta))\)
• The logic of the WLS estimation method is to weigh the correlation elements by the inverse of its sampling covariance matrix.

Fixed-effects TSSEM: Stage 2 analysis (3)

• An example:
  • The sample sizes for the correlations are 200, 500, and 1,000.
  • The asymptotic variance on \(\bar{r}_{21}\) is larger than those on \(\bar{r}_{31}\) and \(\bar{r}_{32}\) because \(\bar{r}_{21}\) is based on smaller sample size (and thus, a larger SE).
  • Because \(V_\mathrm{F}\) is inverted, less weight is given to \(\bar{r}_{21}\) than to \(\bar{r}_{32}\). Therefore, the TSSEM approach takes the precision of the estimates from the stage one analysis into account by using different weights.

## Asymptotic sampling covariance matrix based on the harmonic mean (375):

##           x2_x1     x3_x1     x3_x2
## x2_x1 0.0010923 0.0004864 0.0004864
## x3_x1 0.0004864 0.0010923 0.0004864
## x3_x2 0.0004864 0.0004864 0.0010923

## Asymptotic sampling covariance matrix calculated in TSSEM approach:

##           x2_x1     x3_x1     x3_x2
## x2_x1 0.0020480 0.0005768 0.0004079
## x3_x1 0.0005768 0.0008192 0.0002580
## x3_x2 0.0004079 0.0002580 0.0004096

Fixed-effects TSSEM: Stage 2 analysis (4)

• When a correlation matrix is used with WLS, we have to ensure that the diagonals of the model implied matrix are 1; otherwise, the results are not correct.
• The metaSEM package provides two methods to achieve this:
  • Exclude the error variances in the analysis and calculate them after the analysis (diag.constraints=FALSE);
  • Impose nonlinear constraints to ensure the diagonals are 1 (diag.constraints=TRUE)

Subgroup analysis

• When the assumption of the common correlation matrix is not tenable, we may try two alternative approaches: A subgroup analysis or random-effects model
• Subgroup analysis: try to group studies according to the study characteristics into relatively homogeneous subgroups (Cheung & Chan, 2005; Jak & Cheung, 2018);¹⁶, ¹⁷
  • Pros: Similar to conventional meta-regression by using study characteristics as the moderators;
  • Cons: Only categorical moderators are allowed; few studies in each group.

Random-effects TSSEM

• The fixed-effects TSSEM has been extended to a random-effects model (Cheung, 2014), which is mathematically the same as the random-effects GLS proposed by Becker (1992).
• There are several pros and cons of applying a random-effects model:
  • Pros: Treating differences across studies as random effects;
  • Cons: Difficult to explain study differences by using study characteristics as the moderators.
• Stage 1 analysis:
  • Multivariate random-effects model is used to synthesize correlation matrices;
• Stage 2 analysis:
  • WLS estimation is used to fit structural equation models on the average correlation matrix.

Random-effects TSSEM: Stage 1 analysis (1)

• A multivariate meta-analysis is used to synthesize the correlation matrices:
  • \(\begin{split} \mbox{Level 1: } & r_i = \rho_i + e_i, \\ \mbox{Level 2: } & \rho_i = \rho_\mathrm{R} + u_i, \end{split}\),
  • where \(e_i \sim \mathcal{N}(0, V_i)\) is the known sampling covariance matrix and \(u_i \sim \mathcal{N}(0, T^2)\) is the heterogeneity covariance matrix that has to be estimated.
• The SEM-based multivariate meta-analysis is used to obtain the average correlation matrix \(R_\mathrm{R}\) and its asymptotic sampling covariance matrix \(V_\mathrm{R}\).

Random-effects TSSEM: Stage 1 analysis (2)

• \(I^2\) can be calculated for each correlation coefficients to indicate the degree of heterogeneity.
  
• Notes on the heterogeneity matrix (\(T^2\)) in stage 1 analysis:
  • The no. of effect sizes is usually much larger in MASEM than in conventional multivariate meta-analysis.
  • For example, it is uncommon to have \(p=5\) variables in a MASEM. There will be \(5 \times (5-1)/2=10\) effect sizes and \(10 \times (10+1)/2=55\) elements in \(T^2\).
  • There may not be enough studies to estimate the full \(T^2\).
  • If this happens, we may use a diagonal matrix for \(T^2\) (RE.type="Diag").

Random-effects TSSEM: Stage 2 analysis

• Suppose the proposed correlation structure is \(\rho_\mathrm{R}(\theta)=vechs(P_\mathrm{R}(\theta))\), the discrepancy function is
  • \(F_\mathrm{WLS}(\theta) = (r_\mathrm{R} - \rho_\mathrm{R}(\theta))^T V_\mathrm{R}^{-1} (r_\mathrm{R} - \rho_\mathrm{R}(\theta))\)
  • The model is similar to that for the fixed-effects TSSEM.
• Notes on the heterogeneity matrix (\(T^2\)) in stage 2 analysis:
  • \(T^2\) is not directly involved in the above fit function;
  • As \(V_\mathrm{R}\) is estimated after controlling for the random effects \(T^2\), \(V_\mathrm{R}\) has already been taken the random effects into account;
  • \(V_\mathrm{R}\) is usually larger than \(V_\mathrm{F}\). Thus, the SEs in the random-effects model are generally larger than those in the fixed-effects model.

Correction for unreliability

• The observed correlation coefficients are usually smaller than the actual correlations without measurement error.
• There is some controversy over whether it is necessary to correct for attenuation and the statistical artifacts.
• Rosenthal (1991) criticized the use of a correction for attenuation because the corrected values are not as useful as the uncorrected values in realistic settings.¹⁸
• Other researchers (e.g., Schmidt and Hunter, 2015) argued for correcting for attenuation before combining the correlation matrices.
• MASEM:¹⁹
  • Michel et al. (2011) suggested that substantive model conclusions in the psychological literature are generally unaffected by study artifacts and related statistical corrections.²⁰
  • As their conclusions are based on real examples rather than on computer simulations, further research may address the effects of unreliability in MASEM.

Illustrations with R

• In this section, we are going to illustrate the TSSEM using the metaSEM package.
  • Stage 1 analysis: tssem1()
  • Stage 2 analysis: tssem2()
  • Note: The tssem2() automatically handles whether a fixed- or random-effects model is used in the stage 1 analysis.
• Structural equation models are specified via the reticular action model (RAM) specification (McArdle & McDonald, 1984).²¹
• The metaSEM package also provides a function lavaan2RAM() to convert the lavaan syntax into RAM specification.
• Besides using the likelihood ratio (LR) statistic to assessing the exact fit, we may also use the RMSEA to asses the approximate fit (Browne & Cudeck, 1993):²²
  • If the RMSEA is close to 0.05, the proposed model is considered fitting the data reasonably well;
  • If the RMSEA is larger than 0.10, the proposed model does not fit the data.

A higher-order CFA model for the Big Five model

• Digman (1997) reported 14 correlation matrices among the five-factor model.²³
• Alpha factor (general personality of socialization)
  • agreeableness (A), conscientiousness (C), and emotional stability (ES)
• Beta factor (personal growth)
  • extraversion (E), and intellect (I)

The dataset

• The dataset is available in the metaSEM package as Digman97.
• The same data in Excel format: Digman97.xlsx
• The first sheet (Info) includes the study characteristics whereas the other sheets (Study 1 to Study 14) store the correlation matrices.

Reading the Excel file into R

## Load the library to read XLSX file
library(readxl)

## Read the study characteristics
study <- read_xlsx("Digman97.xlsx", sheet="Info")

head(study)

## # A tibble: 6 x 3
##   Study                               n Cluster     
##   <chr>                           <dbl> <chr>       
## 1 Digman 1 (1994)                   102 Children    
## 2 Digman 2 (1994)                   149 Children    
## 3 Digman 3 (1963c)                  334 Children    
## 4 Digman & Takemoto-Chock (1981b)   162 Children    
## 5 Graziano & Ward (1992)             91 Adolescents 
## 6 Yik & Bond (1993)                 656 Young adults

## Create an empty list to store the correlation matrices
Digman97.data <- list()
  
## Read 1 to 14 correlation matrices
for (i in 1:14) {
  
  ## Read each sheet and convert it into a matrix
  mat <- as.matrix(read_xlsx("Digman97.xlsx", sheet=paste0("Study ", i)))
  
  ## Add the row names
  rownames(mat) <- colnames(mat)
  
  ## Save it into a list
  Digman97.data[[i]] <- mat
}

## Add the names of the studies
names(Digman97.data) <- study$Study

## Show the first few studies
head(Digman97.data)

## $`Digman 1 (1994)`
##        A     C   ES     E    I
## A   1.00  0.62 0.41 -0.48 0.00
## C   0.62  1.00 0.59 -0.10 0.35
## ES  0.41  0.59 1.00  0.27 0.41
## E  -0.48 -0.10 0.27  1.00 0.37
## I   0.00  0.35 0.41  0.37 1.00
## 
## $`Digman 2 (1994)`
##        A    C   ES     E     I
## A   1.00 0.39 0.53 -0.30 -0.05
## C   0.39 1.00 0.59  0.07  0.44
## ES  0.53 0.59 1.00  0.09  0.22
## E  -0.30 0.07 0.09  1.00  0.45
## I  -0.05 0.44 0.22  0.45  1.00
## 
## $`Digman 3 (1963c)`
##       A     C   ES     E    I
## A  1.00  0.65 0.35  0.25 0.14
## C  0.65  1.00 0.37 -0.10 0.33
## ES 0.35  0.37 1.00  0.24 0.41
## E  0.25 -0.10 0.24  1.00 0.41
## I  0.14  0.33 0.41  0.41 1.00
## 
## $`Digman & Takemoto-Chock (1981b)`
##        A     C   ES     E     I
## A   1.00  0.65 0.70 -0.26 -0.03
## C   0.65  1.00 0.71 -0.16  0.24
## ES  0.70  0.71 1.00  0.01  0.11
## E  -0.26 -0.16 0.01  1.00  0.66
## I  -0.03  0.24 0.11  0.66  1.00
## 
## $`Graziano & Ward (1992)`
##       A    C   ES    E    I
## A  1.00 0.64 0.35 0.29 0.22
## C  0.64 1.00 0.27 0.16 0.22
## ES 0.35 0.27 1.00 0.32 0.36
## E  0.29 0.16 0.32 1.00 0.53
## I  0.22 0.22 0.36 0.53 1.00
## 
## $`Yik & Bond (1993)`
##       A    C   ES    E    I
## A  1.00 0.66 0.57 0.35 0.38
## C  0.66 1.00 0.45 0.20 0.31
## ES 0.57 0.45 1.00 0.49 0.31
## E  0.35 0.20 0.49 1.00 0.59
## I  0.38 0.31 0.31 0.59 1.00

## Extract the sample sizes
Digman97.n <- study$n
Digman97.n

##  [1]  102  149  334  162   91  656   70   70  277  227 1000  227   91 1040

## Extract the cluster
Digman97.cluster <- study$Cluster
Digman97.cluster

##  [1] "Children"      "Children"      "Children"      "Children"      "Adolescents"   "Young adults"  "Young adults"  "Young adults"  "Mature adults"
## [10] "Mature adults" "Mature adults" "Mature adults" "Mature adults" "Mature adults"

Proposed model

• RAM model: the observed \((o)\) and latent variables \((l)\) are combined:
  • A matrix: \((o+l) \times (o+l))\) regression paths and factor loadings
  • S matrix: \((o+l) \times (o+l))\) variances, and covariances
  • F matrix: \(o \times (o+l))\) filter matrix to select the observed variables
• Features of the model specification
  • Numbers: fixed values, e.g., 0 and 1 means that the values are fixed at 0 and 1, respectively.
  • "0.3*Alpha_A": the factor loading from Alpha to A is free with a starting value of 0.3.
  • If the labels are the same, these parameters are constrained to be equal.
• Important: The columns and rows represent the independent and dependent variables in the A matrix!
• Most users find it easier to specify the models in lavaan syntax.

model <- "## Factor loadings
          ## Alpha is measured by A, C, and ES
          Alpha =~ A + C + ES
          ## Beta is measured by E and I
          Beta =~ E + I
          ## Factor correlation between Alpha and Beta
          Alpha ~~ Beta"

## Display the model
plot(model, color="yellow")

## Convert the lavaan syntax into a RAM model as the metaSEM only knows the RAM model
## It is important to ensure that the variables are arranged in A, C, ES, E, and I.
RAM <- lavaan2RAM(model, obs.variables=c("A","C","ES","E","I"),
                  A.notation="on", S.notation="with")
RAM

## $A
##       A   C   ES  E   I   Alpha         Beta       
## A     "0" "0" "0" "0" "0" "0*AonAlpha"  "0"        
## C     "0" "0" "0" "0" "0" "0*ConAlpha"  "0"        
## ES    "0" "0" "0" "0" "0" "0*ESonAlpha" "0"        
## E     "0" "0" "0" "0" "0" "0"           "0*EonBeta"
## I     "0" "0" "0" "0" "0" "0"           "0*IonBeta"
## Alpha "0" "0" "0" "0" "0" "0"           "0"        
## Beta  "0" "0" "0" "0" "0" "0"           "0"        
## 
## $S
##       A          C          ES           E          I          Alpha             Beta             
## A     "0*AwithA" "0"        "0"          "0"        "0"        "0"               "0"              
## C     "0"        "0*CwithC" "0"          "0"        "0"        "0"               "0"              
## ES    "0"        "0"        "0*ESwithES" "0"        "0"        "0"               "0"              
## E     "0"        "0"        "0"          "0*EwithE" "0"        "0"               "0"              
## I     "0"        "0"        "0"          "0"        "0*IwithI" "0"               "0"              
## Alpha "0"        "0"        "0"          "0"        "0"        "1"               "0*AlphawithBeta"
## Beta  "0"        "0"        "0"          "0"        "0"        "0*AlphawithBeta" "1"              
## 
## $F
##    A C ES E I Alpha Beta
## A  1 0  0 0 0     0    0
## C  0 1  0 0 0     0    0
## ES 0 0  1 0 0     0    0
## E  0 0  0 1 0     0    0
## I  0 0  0 0 1     0    0
## 
## $M
##   A C ES E I Alpha Beta
## 1 0 0  0 0 0     0    0

Fixed-effects TSSEM: Stage 1 Analysis

• The test statistic on testing the homogeneity test of the correlation matrix is \(\chi^2(df=130, N=4496)=1505.44, p < 0.001\).
• Both the RMSEA and SRMR are very large indicating that the proposed model (homogeneity of correlation matrices) did not fit the data well.
• The estimates of the \(S[i,j]\) represent the element of the pooled correlation (or covariance) matrix.

## method="FEM": fixed-effects TSSEM
fixed1 <- tssem1(Digman97.data, Digman97.n, method="FEM")

## summary of the findings
summary(fixed1)

## 
## Call:
## tssem1FEM(Cov = Cov, n = n, cor.analysis = cor.analysis, model.name = model.name, 
##     cluster = cluster, suppressWarnings = suppressWarnings, silent = silent, 
##     run = run)
## 
## Coefficients:
##        Estimate Std.Error z value  Pr(>|z|)    
## S[1,2] 0.363278  0.013368 27.1760 < 2.2e-16 ***
## S[1,3] 0.390375  0.012880 30.3076 < 2.2e-16 ***
## S[1,4] 0.103572  0.015048  6.8830 5.862e-12 ***
## S[1,5] 0.092286  0.015047  6.1331 8.621e-10 ***
## S[2,3] 0.416070  0.012519 33.2344 < 2.2e-16 ***
## S[2,4] 0.135148  0.014776  9.1464 < 2.2e-16 ***
## S[2,5] 0.141431  0.014866  9.5135 < 2.2e-16 ***
## S[3,4] 0.244459  0.014153 17.2724 < 2.2e-16 ***
## S[3,5] 0.138339  0.014834  9.3259 < 2.2e-16 ***
## S[4,5] 0.424566  0.012375 34.3070 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                      Value
## Sample size                      4496.0000
## Chi-square of target model       1505.4443
## DF of target model                130.0000
## p value of target model             0.0000
## Chi-square of independence model 4471.4242
## DF of independence model          140.0000
## RMSEA                               0.1815
## RMSEA lower 95% CI                  0.1736
## RMSEA upper 95% CI                  0.1901
## SRMR                                0.1621
## TLI                                 0.6580
## CFI                                 0.6824
## AIC                              1245.4443
## BIC                               412.0217
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)

## extract coefficients
coef(fixed1)

##             A         C        ES         E          I
## A  1.00000000 0.3632782 0.3903748 0.1035716 0.09228557
## C  0.36327824 1.0000000 0.4160695 0.1351477 0.14143058
## ES 0.39037483 0.4160695 1.0000000 0.2444593 0.13833895
## E  0.10357155 0.1351477 0.2444593 1.0000000 0.42456626
## I  0.09228557 0.1414306 0.1383390 0.4245663 1.00000000

Fixed-effects TSSEM: Stage 2 Analysis (1)

• Although the random-effects model is preferred, we illustrate how to fit the stage 2 analysis based on the fixed-effects model.
• Important notes:
  • The SEs of the parameter estimates are likely underestimated when the correlation matrices are not homogeneous.
  • Readers should be cautious when interpreting the results of the stage 2 analysis when the homogeneity of the correlation matrices is rejected in the stage 1 analysis.

fixed2 <- tssem2(fixed1, Amatrix=RAM$A, Smatrix=RAM$S, Fmatrix=RAM$F)
summary(fixed2)

## 
## Call:
## wls(Cov = coef.tssem1FEM(tssem1.obj), aCov = vcov.tssem1FEM(tssem1.obj), 
##     n = sum(tssem1.obj$n), Amatrix = Amatrix, Smatrix = Smatrix, 
##     Fmatrix = Fmatrix, diag.constraints = diag.constraints, cor.analysis = tssem1.obj$cor.analysis, 
##     intervals.type = intervals.type, mx.algebras = mx.algebras, 
##     model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##               Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
## AonAlpha      0.562800  0.015380 0.532655 0.592944  36.593 < 2.2e-16 ***
## ConAlpha      0.605217  0.015324 0.575183 0.635251  39.495 < 2.2e-16 ***
## EonBeta       0.781413  0.034244 0.714297 0.848529  22.819 < 2.2e-16 ***
## ESonAlpha     0.719201  0.015685 0.688458 0.749943  45.852 < 2.2e-16 ***
## IonBeta       0.551374  0.026031 0.500354 0.602394  21.181 < 2.2e-16 ***
## AlphawithBeta 0.362678  0.022384 0.318806 0.406550  16.203 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                                Value
## Sample size                                4496.0000
## Chi-square of target model                   65.4526
## DF of target model                            4.0000
## p value of target model                       0.0000
## Number of constraints imposed on "Smatrix"    0.0000
## DF manually adjusted                          0.0000
## Chi-square of independence model           3112.7380
## DF of independence model                     10.0000
## RMSEA                                         0.0585
## RMSEA lower 95% CI                            0.0465
## RMSEA upper 95% CI                            0.0713
## SRMR                                          0.0284
## TLI                                           0.9505
## CFI                                           0.9802
## AIC                                          57.4526
## BIC                                          31.8088
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)

Fixed-effects TSSEM: Stage 2 Analysis (2)

• We may plot the figure using the semPlot package (Epskamp, 2014).²⁴

plot(fixed2, color="green")

Fixed-effects TSSEM with cluster: Stage 1 Analysis

• As the assumption of the homogeneity of the correlation matrices has not been met, we group the studies into several clusters based on the study characteristics.
• For ease of illustration, we recode the data into younger and older participants.
• From the results, splitting the data into younger and older participants does not explain the heterogeneity of the data well.

# Display the original study characteristic
table(Digman97.cluster)     

## Digman97.cluster
##   Adolescents      Children Mature adults  Young adults 
##             1             4             6             3

## Younger participants: "Children" and "Adolescents"
## Older participants: "Mature adults"
sample <- ifelse(Digman97.cluster %in% c("Children", "Adolescents"), 
                 yes="Younger participants", no="Older participants")

table(sample)

## sample
##   Older participants Younger participants 
##                    9                    5

## cluster: variable for the analysis with cluster
fixed1.cluster <- tssem1(Digman97.data, Digman97.n, method="FEM", cluster=sample)

summary(fixed1.cluster)

## $`Older participants`
## 
## Call:
## tssem1FEM(Cov = data.cluster[[i]], n = n.cluster[[i]], cor.analysis = cor.analysis, 
##     model.name = model.name, suppressWarnings = suppressWarnings)
## 
## Coefficients:
##        Estimate Std.Error z value  Pr(>|z|)    
## S[1,2] 0.297471  0.015436 19.2716 < 2.2e-16 ***
## S[1,3] 0.370248  0.014532 25.4774 < 2.2e-16 ***
## S[1,4] 0.137694  0.016403  8.3944 < 2.2e-16 ***
## S[1,5] 0.098061  0.016724  5.8637 4.528e-09 ***
## S[2,3] 0.393692  0.014146 27.8306 < 2.2e-16 ***
## S[2,4] 0.183045  0.016055 11.4009 < 2.2e-16 ***
## S[2,5] 0.092774  0.016643  5.5743 2.485e-08 ***
## S[3,4] 0.260753  0.015554 16.7645 < 2.2e-16 ***
## S[3,5] 0.096141  0.016573  5.8009 6.596e-09 ***
## S[4,5] 0.411756  0.013900 29.6225 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                      Value
## Sample size                      3658.0000
## Chi-square of target model        825.9826
## DF of target model                 80.0000
## p value of target model             0.0000
## Chi-square of independence model 3000.9731
## DF of independence model           90.0000
## RMSEA                               0.1515
## RMSEA lower 95% CI                  0.1424
## RMSEA upper 95% CI                  0.1611
## SRMR                                0.1459
## TLI                                 0.7117
## CFI                                 0.7437
## AIC                               665.9826
## BIC                               169.6088
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## 
## $`Younger participants`
## 
## Call:
## tssem1FEM(Cov = data.cluster[[i]], n = n.cluster[[i]], cor.analysis = cor.analysis, 
##     model.name = model.name, suppressWarnings = suppressWarnings)
## 
## Coefficients:
##         Estimate Std.Error z value  Pr(>|z|)    
## S[1,2]  0.604330  0.022125 27.3141 < 2.2e-16 ***
## S[1,3]  0.465536  0.027493 16.9327 < 2.2e-16 ***
## S[1,4] -0.031381  0.035940 -0.8732   0.38258    
## S[1,5]  0.061508  0.034547  1.7804   0.07500 .  
## S[2,3]  0.501432  0.026348 19.0311 < 2.2e-16 ***
## S[2,4] -0.060678  0.034557 -1.7559   0.07911 .  
## S[2,5]  0.320987  0.031064 10.3330 < 2.2e-16 ***
## S[3,4]  0.175437  0.033675  5.2097 1.891e-07 ***
## S[3,5]  0.305149  0.031586  9.6609 < 2.2e-16 ***
## S[4,5]  0.478640  0.026883 17.8045 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                      Value
## Sample size                       838.0000
## Chi-square of target model        346.2810
## DF of target model                 40.0000
## p value of target model             0.0000
## Chi-square of independence model 1470.4511
## DF of independence model           50.0000
## RMSEA                               0.2139
## RMSEA lower 95% CI                  0.1939
## RMSEA upper 95% CI                  0.2355
## SRMR                                0.1411
## TLI                                 0.7305
## CFI                                 0.7844
## AIC                               266.2810
## BIC                                77.0402
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)

Fixed-effects TSSEM with cluster: Stage 2 Analysis (1)

• As an illustration, we show how to fit the stage 2 analysis.
• We should not trust the results because the estimates, e.g., Beta_I=3.29, for the younger participants are problematic.

fixed2.cluster <- tssem2(fixed1.cluster, Amatrix=RAM$A, Smatrix=RAM$S, Fmatrix=RAM$F)

summary(fixed2.cluster)

## $`Older participants`
## 
## Call:
## wls(Cov = coef.tssem1FEM(tssem1.obj), aCov = vcov.tssem1FEM(tssem1.obj), 
##     n = sum(tssem1.obj$n), Amatrix = Amatrix, Smatrix = Smatrix, 
##     Fmatrix = Fmatrix, diag.constraints = diag.constraints, cor.analysis = tssem1.obj$cor.analysis, 
##     intervals.type = intervals.type, mx.algebras = mx.algebras, 
##     model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##               Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
## AonAlpha      0.512656  0.018206 0.476973 0.548340  28.158 < 2.2e-16 ***
## ConAlpha      0.549967  0.017945 0.514795 0.585140  30.647 < 2.2e-16 ***
## EonBeta       0.967136  0.058751 0.851986 1.082287  16.462 < 2.2e-16 ***
## ESonAlpha     0.732174  0.018929 0.695074 0.769273  38.680 < 2.2e-16 ***
## IonBeta       0.430649  0.029634 0.372568 0.488730  14.532 < 2.2e-16 ***
## AlphawithBeta 0.349236  0.028118 0.294125 0.404346  12.420 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                                Value
## Sample size                                3658.0000
## Chi-square of target model                   21.9954
## DF of target model                            4.0000
## p value of target model                       0.0002
## Number of constraints imposed on "Smatrix"    0.0000
## DF manually adjusted                          0.0000
## Chi-square of independence model           2273.3236
## DF of independence model                     10.0000
## RMSEA                                         0.0351
## RMSEA lower 95% CI                            0.0217
## RMSEA upper 95% CI                            0.0500
## SRMR                                          0.0160
## TLI                                           0.9801
## CFI                                           0.9920
## AIC                                          13.9954
## BIC                                         -10.8233
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
## 
## $`Younger participants`
## 
## Call:
## wls(Cov = coef.tssem1FEM(tssem1.obj), aCov = vcov.tssem1FEM(tssem1.obj), 
##     n = sum(tssem1.obj$n), Amatrix = Amatrix, Smatrix = Smatrix, 
##     Fmatrix = Fmatrix, diag.constraints = diag.constraints, cor.analysis = tssem1.obj$cor.analysis, 
##     intervals.type = intervals.type, mx.algebras = mx.algebras, 
##     model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##                Estimate Std.Error    lbound    ubound z value Pr(>|z|)    
## AonAlpha       0.747647  0.023880  0.700842  0.794451 31.3081   <2e-16 ***
## ConAlpha       0.911705  0.019864  0.872772  0.950638 45.8968   <2e-16 ***
## EonBeta        0.152561  0.159140 -0.159347  0.464469  0.9587   0.3377    
## ESonAlpha      0.677434  0.025864  0.626742  0.728126 26.1924   <2e-16 ***
## IonBeta        3.283875  3.363587 -3.308633  9.876384  0.9763   0.3289    
## AlphawithBeta  0.117255  0.125430 -0.128584  0.363094  0.9348   0.3499    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                                Value
## Sample size                                 838.0000
## Chi-square of target model                  145.6167
## DF of target model                            4.0000
## p value of target model                       0.0000
## Number of constraints imposed on "Smatrix"    0.0000
## DF manually adjusted                          0.0000
## Chi-square of independence model           2480.2352
## DF of independence model                     10.0000
## RMSEA                                         0.2057
## RMSEA lower 95% CI                            0.1778
## RMSEA upper 95% CI                            0.2350
## SRMR                                          0.1051
## TLI                                           0.8567
## CFI                                           0.9427
## AIC                                         137.6167
## BIC                                         118.6926
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)

Fixed-effects TSSEM with cluster: Stage 2 Analysis (2)

• The multiple-group analysis may also be plotted.

## Setup two plots
layout(t(1:2))
invisible(lapply(fixed2.cluster, plot))

• We may also plot the figures with more details

## Setup two plots
layout(t(1:2))

## Plot the first group
plot(fixed2.cluster[[1]])
title("Younger participants")

## Plot the second group
plot(fixed2.cluster[[2]])
title("Older participants")

Random-effects TSSEM: Stage 1 Analysis

• There is a total of ten correlation coefficients in the analysis with 55 elements in the variance component matrix of the random effects.

• We may specify RE.type="Diag" with only 10 variances in the diagonal variance component matrix.

• The minimum and the maximum \(I^2\) are 0.7669 and 0.9326, respectively.

• The random-effects model is more appropriate.

## method="REM": Random-effects model
random1 <- tssem1(Digman97.data, Digman97.n, method="REM", RE.type="Diag")

summary(random1)

## 
## Call:
## meta(y = ES, v = acovR, RE.constraints = Diag(paste0(RE.startvalues, 
##     "*Tau2_", 1:no.es, "_", 1:no.es)), RE.lbound = RE.lbound, 
##     I2 = I2, model.name = model.name, suppressWarnings = TRUE, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##                Estimate   Std.Error      lbound      ubound z value  Pr(>|z|)    
## Intercept1   0.38971902  0.05429255  0.28330757  0.49613047  7.1781 7.068e-13 ***
## Intercept2   0.43265876  0.04145135  0.35141562  0.51390190 10.4377 < 2.2e-16 ***
## Intercept3   0.04945623  0.06071078 -0.06953471  0.16844718  0.8146   0.41529    
## Intercept4   0.09603703  0.04478710  0.00825593  0.18381813  2.1443   0.03201 *  
## Intercept5   0.42724236  0.03911620  0.35057603  0.50390870 10.9224 < 2.2e-16 ***
## Intercept6   0.11929312  0.04106203  0.03881303  0.19977322  2.9052   0.00367 ** 
## Intercept7   0.19292424  0.04757961  0.09966991  0.28617856  4.0548 5.018e-05 ***
## Intercept8   0.22690158  0.03240893  0.16338125  0.29042191  7.0012 2.538e-12 ***
## Intercept9   0.18105561  0.04258855  0.09758360  0.26452763  4.2513 2.126e-05 ***
## Intercept10  0.43614966  0.03205960  0.37331400  0.49898532 13.6043 < 2.2e-16 ***
## Tau2_1_1     0.03648370  0.01513054  0.00682839  0.06613902  2.4113   0.01590 *  
## Tau2_2_2     0.01963096  0.00859788  0.00277942  0.03648250  2.2832   0.02242 *  
## Tau2_3_3     0.04571436  0.01952284  0.00745029  0.08397842  2.3416   0.01920 *  
## Tau2_4_4     0.02236119  0.00995082  0.00285794  0.04186445  2.2472   0.02463 *  
## Tau2_5_5     0.01729071  0.00796404  0.00168149  0.03289994  2.1711   0.02992 *  
## Tau2_6_6     0.01815481  0.00895896  0.00059557  0.03571405  2.0264   0.04272 *  
## Tau2_7_7     0.02604881  0.01130265  0.00389601  0.04820160  2.3047   0.02119 *  
## Tau2_8_8     0.00988761  0.00513713 -0.00018098  0.01995619  1.9247   0.05426 .  
## Tau2_9_9     0.01988243  0.00895052  0.00233973  0.03742513  2.2214   0.02633 *  
## Tau2_10_10   0.01049221  0.00501578  0.00066147  0.02032295  2.0918   0.03645 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 1220.332
## Degrees of freedom of the Q statistic: 130
## P value of the Q statistic: 0
## 
## Heterogeneity indices (based on the estimated Tau2):
##                               Estimate
## Intercept1: I2 (Q statistic)    0.9326
## Intercept2: I2 (Q statistic)    0.8872
## Intercept3: I2 (Q statistic)    0.9325
## Intercept4: I2 (Q statistic)    0.8703
## Intercept5: I2 (Q statistic)    0.8797
## Intercept6: I2 (Q statistic)    0.8480
## Intercept7: I2 (Q statistic)    0.8887
## Intercept8: I2 (Q statistic)    0.7669
## Intercept9: I2 (Q statistic)    0.8590
## Intercept10: I2 (Q statistic)   0.8193
## 
## Number of studies (or clusters): 14
## Number of observed statistics: 140
## Number of estimated parameters: 20
## Degrees of freedom: 120
## -2 log likelihood: -112.4196 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)

## Extract the fixed-effects estimates
(est_fixed <- coef(random1, select="fixed"))

##  Intercept1  Intercept2  Intercept3  Intercept4  Intercept5  Intercept6  Intercept7  Intercept8  Intercept9 Intercept10 
##  0.38971902  0.43265876  0.04945623  0.09603703  0.42724236  0.11929312  0.19292424  0.22690158  0.18105561  0.43614966

## Convert the estimated vector to a symmetrical matrix
## where the diagonals are fixed at 1 (for a correlation matrix)
vec2symMat(est_fixed, diag=FALSE)

##            [,1]      [,2]      [,3]       [,4]       [,5]
## [1,] 1.00000000 0.3897190 0.4326588 0.04945623 0.09603703
## [2,] 0.38971902 1.0000000 0.4272424 0.11929312 0.19292424
## [3,] 0.43265876 0.4272424 1.0000000 0.22690158 0.18105561
## [4,] 0.04945623 0.1192931 0.2269016 1.00000000 0.43614966
## [5,] 0.09603703 0.1929242 0.1810556 0.43614966 1.00000000

Random-effects TSSEM: Stage 2 Analysis

• The stage 2 analysis is similar to that for the fixed-effects model.

• We may trust the results based on the random-effects model.
• The proposed model fits the data well with \(\chi^2(df=4, N=4496)=7.82, p=0.0984\), CFI=\(0.9922\), RMSEA=\(0.0146\), and SRMR=\(0.0436\).

• The estimated correlation between the two factors (labeled cor in the printout) is \(\hat{S}[7,6]=0.3777\).

random2 <- tssem2(random1, Amatrix=RAM$A, Smatrix=RAM$S, Fmatrix=RAM$F)
summary(random2)

## 
## Call:
## wls(Cov = pooledS, aCov = aCov, n = tssem1.obj$total.n, Amatrix = Amatrix, 
##     Smatrix = Smatrix, Fmatrix = Fmatrix, diag.constraints = diag.constraints, 
##     cor.analysis = cor.analysis, intervals.type = intervals.type, 
##     mx.algebras = mx.algebras, model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##               Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
## AonAlpha      0.569435  0.052425 0.466684 0.672187 10.8619 < 2.2e-16 ***
## ConAlpha      0.590630  0.052649 0.487439 0.693821 11.2182 < 2.2e-16 ***
## EonBeta       0.679955  0.075723 0.531541 0.828370  8.9795 < 2.2e-16 ***
## ESonAlpha     0.760455  0.061963 0.639009 0.881901 12.2726 < 2.2e-16 ***
## IonBeta       0.641842  0.072459 0.499825 0.783859  8.8580 < 2.2e-16 ***
## AlphawithBeta 0.377691  0.047402 0.284785 0.470596  7.9679 1.554e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                                Value
## Sample size                                4496.0000
## Chi-square of target model                    7.8204
## DF of target model                            4.0000
## p value of target model                       0.0984
## Number of constraints imposed on "Smatrix"    0.0000
## DF manually adjusted                          0.0000
## Chi-square of independence model            501.6767
## DF of independence model                     10.0000
## RMSEA                                         0.0146
## RMSEA lower 95% CI                            0.0000
## RMSEA upper 95% CI                            0.0297
## SRMR                                          0.0436
## TLI                                           0.9806
## CFI                                           0.9922
## AIC                                          -0.1796
## BIC                                         -25.8234
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)

Plotting the results

## Plot the parameter estimates
plot(random2, color="green")

A regression model on SAT (Math) (1)

• The dataset is available in the metaSEM package as Becker94.
• The data file in Excel format: Becker94.xlsx

## Read the first sheet in the Excel file
my.xlsx <- read_excel("Becker94.xlsx")
my.xlsx

## # A tibble: 10 x 6
##    Study                    r_SAT_Math_Spatial r_SAT_Math_SAT_Verbal Spatial_SAT_Verbal     n gender 
##    <chr>                                 <dbl>                 <dbl>              <dbl> <dbl> <chr>  
##  1 Becker (1978) Females                 0.47                  -0.21              -0.15    74 Females
##  2 Becker (1978) Males                   0.28                   0.19               0.18   153 Males  
##  3 Berry (1957) Females                  0.48                   0.41               0.26    48 Females
##  4 Berry (1957) Males                    0.37                   0.4                0.27    55 Males  
##  5 Rosenberg (1981) Females              0.42                   0.48               0.23    51 Females
##  6 Rosenberg (1981) Males                0.41                   0.74               0.44    18 Males  
##  7 Weiner A (1984) Females               0.26                   0.72               0.36    27 Females
##  8 Weiner A (1984) Males                 0.32                   0.52               0.1     43 Males  
##  9 Weiner B (1984) Females               0.580                  0.64               0.4     35 Females
## 10 Weiner B (1984) Males                 0.34                   0.28              -0.03    34 Males

## Split the data by rows
my.R2 <- split(my.xlsx[, 2:4], seq(nrow(my.xlsx)))
head(my.R2)

## $`1`
## # A tibble: 1 x 3
##   r_SAT_Math_Spatial r_SAT_Math_SAT_Verbal Spatial_SAT_Verbal
##                <dbl>                 <dbl>              <dbl>
## 1               0.47                 -0.21              -0.15
## 
## $`2`
## # A tibble: 1 x 3
##   r_SAT_Math_Spatial r_SAT_Math_SAT_Verbal Spatial_SAT_Verbal
##                <dbl>                 <dbl>              <dbl>
## 1               0.28                  0.19               0.18
## 
## $`3`
## # A tibble: 1 x 3
##   r_SAT_Math_Spatial r_SAT_Math_SAT_Verbal Spatial_SAT_Verbal
##                <dbl>                 <dbl>              <dbl>
## 1               0.48                  0.41               0.26
## 
## $`4`
## # A tibble: 1 x 3
##   r_SAT_Math_Spatial r_SAT_Math_SAT_Verbal Spatial_SAT_Verbal
##                <dbl>                 <dbl>              <dbl>
## 1               0.37                   0.4               0.27
## 
## $`5`
## # A tibble: 1 x 3
##   r_SAT_Math_Spatial r_SAT_Math_SAT_Verbal Spatial_SAT_Verbal
##                <dbl>                 <dbl>              <dbl>
## 1               0.42                  0.48               0.23
## 
## $`6`
## # A tibble: 1 x 3
##   r_SAT_Math_Spatial r_SAT_Math_SAT_Verbal Spatial_SAT_Verbal
##                <dbl>                 <dbl>              <dbl>
## 1               0.41                  0.74               0.44

## Convert the row correlations into correlation matrices
my.R2 <- lapply(my.R2, function(x) vec2symMat(unlist(x), diag=FALSE))
head(my.R2)

## $`1`
##       [,1]  [,2]  [,3]
## [1,]  1.00  0.47 -0.21
## [2,]  0.47  1.00 -0.15
## [3,] -0.21 -0.15  1.00
## 
## $`2`
##      [,1] [,2] [,3]
## [1,] 1.00 0.28 0.19
## [2,] 0.28 1.00 0.18
## [3,] 0.19 0.18 1.00
## 
## $`3`
##      [,1] [,2] [,3]
## [1,] 1.00 0.48 0.41
## [2,] 0.48 1.00 0.26
## [3,] 0.41 0.26 1.00
## 
## $`4`
##      [,1] [,2] [,3]
## [1,] 1.00 0.37 0.40
## [2,] 0.37 1.00 0.27
## [3,] 0.40 0.27 1.00
## 
## $`5`
##      [,1] [,2] [,3]
## [1,] 1.00 0.42 0.48
## [2,] 0.42 1.00 0.23
## [3,] 0.48 0.23 1.00
## 
## $`6`
##      [,1] [,2] [,3]
## [1,] 1.00 0.41 0.74
## [2,] 0.41 1.00 0.44
## [3,] 0.74 0.44 1.00

## Add the labels
var.names <- c("SAT_Math", "Spatial", "SAT_Verbal")
my.R2 <- lapply( my.R2, function(x) { dimnames(x) <- list(var.names, var.names); x} )

## Add the study name
names(my.R2) <- my.xlsx$Study
head(my.R2)

## $`Becker (1978) Females`
##            SAT_Math Spatial SAT_Verbal
## SAT_Math       1.00    0.47      -0.21
## Spatial        0.47    1.00      -0.15
## SAT_Verbal    -0.21   -0.15       1.00
## 
## $`Becker (1978) Males`
##            SAT_Math Spatial SAT_Verbal
## SAT_Math       1.00    0.28       0.19
## Spatial        0.28    1.00       0.18
## SAT_Verbal     0.19    0.18       1.00
## 
## $`Berry (1957) Females`
##            SAT_Math Spatial SAT_Verbal
## SAT_Math       1.00    0.48       0.41
## Spatial        0.48    1.00       0.26
## SAT_Verbal     0.41    0.26       1.00
## 
## $`Berry (1957) Males`
##            SAT_Math Spatial SAT_Verbal
## SAT_Math       1.00    0.37       0.40
## Spatial        0.37    1.00       0.27
## SAT_Verbal     0.40    0.27       1.00
## 
## $`Rosenberg (1981) Females`
##            SAT_Math Spatial SAT_Verbal
## SAT_Math       1.00    0.42       0.48
## Spatial        0.42    1.00       0.23
## SAT_Verbal     0.48    0.23       1.00
## 
## $`Rosenberg (1981) Males`
##            SAT_Math Spatial SAT_Verbal
## SAT_Math       1.00    0.41       0.74
## Spatial        0.41    1.00       0.44
## SAT_Verbal     0.74    0.44       1.00

A regression model on SAT (Math) (2)

• Becker and Schram (1994) reported 10 independent samples (five for males and five for females) from five studies with a total sample size of 538.²⁵
• Their study includes correlations among Math, Verbal, and Spatial.
• 

## Regression model
model <- "SAT_Math ~ Spatial + SAT_Verbal
          ## Correlation between Spatial and Verbal
          Spatial ~~ SAT_Verbal
          ## Fix the variances of the independent variables at 1.0
          Spatial ~~ 1*Spatial
          SAT_Verbal ~~ 1*SAT_Verbal"

## Plot the model
plot(model, color="yellow")

## Convert the equation into RAM 
RAM <- lavaan2RAM(model, obs.variables=c("SAT_Math", "Spatial", "SAT_Verbal"))
RAM

## $A
##            SAT_Math Spatial               SAT_Verbal              
## SAT_Math   "0"      "0*SAT_MathONSpatial" "0*SAT_MathONSAT_Verbal"
## Spatial    "0"      "0"                   "0"                     
## SAT_Verbal "0"      "0"                   "0"                     
## 
## $S
##            SAT_Math                 Spatial                   SAT_Verbal               
## SAT_Math   "0*SAT_MathWITHSAT_Math" "0"                       "0"                      
## Spatial    "0"                      "1"                       "0*SpatialWITHSAT_Verbal"
## SAT_Verbal "0"                      "0*SpatialWITHSAT_Verbal" "1"                      
## 
## $F
##            SAT_Math Spatial SAT_Verbal
## SAT_Math          1       0          0
## Spatial           0       1          0
## SAT_Verbal        0       0          1
## 
## $M
##   SAT_Math Spatial SAT_Verbal
## 1        0       0          0

Fixed-effects TSSEM: Stage 1 Analysis

• The test statistic on testing the homogeneity test of the correlation matrix is \(\chi^2(df=27, N=538)=63.66, p < 0.001\).
• Both the RMSEA and SRMR are very large indicating that the proposed model (homogeneity of correlation matrices) did not fit the data well.
• Since the fixed-effects model is not appropriate, we do not fit the stage 2 model.

## method="FEM": fixed-effects TSSEM
fixed1 <- tssem1(Becker94$data, Becker94$n, method="FEM")

## summary of the findings
summary(fixed1)

## 
## Call:
## tssem1FEM(Cov = Cov, n = n, cor.analysis = cor.analysis, model.name = model.name, 
##     cluster = cluster, suppressWarnings = suppressWarnings, silent = silent, 
##     run = run)
## 
## Coefficients:
##        Estimate Std.Error z value  Pr(>|z|)    
## S[1,2] 0.379961  0.037123 10.2351 < 2.2e-16 ***
## S[1,3] 0.334562  0.039947  8.3751 < 2.2e-16 ***
## S[2,3] 0.176461  0.042334  4.1683 3.069e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                      Value
## Sample size                       538.0000
## Chi-square of target model         63.6553
## DF of target model                 27.0000
## p value of target model             0.0001
## Chi-square of independence model  207.7894
## DF of independence model           30.0000
## RMSEA                               0.1590
## RMSEA lower 95% CI                  0.1096
## RMSEA upper 95% CI                  0.2117
## SRMR                                0.1586
## TLI                                 0.7709
## CFI                                 0.7938
## AIC                                 9.6553
## BIC                              -106.1169
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)

Random-effects TSSEM: Stage 1 Analysis

• There is a total of 3 correlation coefficients in the analysis with 6 elements in the variance component matrix of the random effects.

• We may specify RE.type="Diag" since there are only ten studies.

• The minimum and the maximum \(I^2\) are 0.0337 (cor between SAT (Math) and Spatial ability) and 0.7224 (cor between SAT (Math) and SAT (Verbal)), respectively.

• The random-effects model is more appropriate.

## method="REM": Random-effects model
random1 <- tssem1(Becker94$data, Becker94$n, method="REM", RE.type="Diag")
summary(random1)

## 
## Call:
## meta(y = ES, v = acovR, RE.constraints = Diag(paste0(RE.startvalues, 
##     "*Tau2_", 1:no.es, "_", 1:no.es)), RE.lbound = RE.lbound, 
##     I2 = I2, model.name = model.name, suppressWarnings = TRUE, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##              Estimate  Std.Error     lbound     ubound z value  Pr(>|z|)    
## Intercept1  0.3777491  0.0395030  0.3003246  0.4551736  9.5625 < 2.2e-16 ***
## Intercept2  0.3807843  0.0784956  0.2269357  0.5346329  4.8510 1.228e-06 ***
## Intercept3  0.1704927  0.0513545  0.0698398  0.2711457  3.3199 0.0009004 ***
## Tau2_1_1    0.0005038  0.0042009 -0.0077298  0.0087374  0.1199 0.9045407    
## Tau2_2_2    0.0416264  0.0257388 -0.0088206  0.0920734  1.6173 0.1058209    
## Tau2_3_3    0.0067540  0.0102792 -0.0133928  0.0269008  0.6571 0.5111466    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 61.02636
## Degrees of freedom of the Q statistic: 27
## P value of the Q statistic: 0.0001932118
## 
## Heterogeneity indices (based on the estimated Tau2):
##                              Estimate
## Intercept1: I2 (Q statistic)   0.0337
## Intercept2: I2 (Q statistic)   0.7224
## Intercept3: I2 (Q statistic)   0.2676
## 
## Number of studies (or clusters): 10
## Number of observed statistics: 30
## Number of estimated parameters: 6
## Degrees of freedom: 24
## -2 log likelihood: -22.61046 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)

## Extract the fixed-effects estimates
(est_fixed <- coef(random1, select="fixed"))

## Intercept1 Intercept2 Intercept3 
##  0.3777491  0.3807843  0.1704927

## Convert the estimated vector to a symmetrical matrix
## where the diagonals are fixed at 1 (for a correlation matrix)
vec2symMat(est_fixed, diag=FALSE)

##           [,1]      [,2]      [,3]
## [1,] 1.0000000 0.3777491 0.3807843
## [2,] 0.3777491 1.0000000 0.1704927
## [3,] 0.3807843 0.1704927 1.0000000

Random-effects TSSEM: Stage 2 Analysis

• The proposed model is saturated. Thus, the LR statistic is 0.
• The estimated error variance on SAT (Math) is 0.7542. Thus, the \(R^2\) is 0.2458.

random2 <- tssem2(random1, Amatrix=RAM$A, Smatrix=RAM$S, Fmatrix=RAM$F)
summary(random2)

## 
## Call:
## wls(Cov = pooledS, aCov = aCov, n = tssem1.obj$total.n, Amatrix = Amatrix, 
##     Smatrix = Smatrix, Fmatrix = Fmatrix, diag.constraints = diag.constraints, 
##     cor.analysis = cor.analysis, intervals.type = intervals.type, 
##     mx.algebras = mx.algebras, model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##                       Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
## SAT_MathONSAT_Verbal  0.325853  0.080000 0.169056 0.482649  4.0732 4.638e-05 ***
## SAT_MathONSpatial     0.322194  0.042415 0.239063 0.405325  7.5963 3.042e-14 ***
## SpatialWITHSAT_Verbal 0.170493  0.051355 0.069840 0.271146  3.3199 0.0009004 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                             Value
## Sample size                                538.00
## Chi-square of target model                   0.00
## DF of target model                           0.00
## p value of target model                      0.00
## Number of constraints imposed on "Smatrix"   0.00
## DF manually adjusted                         0.00
## Chi-square of independence model           110.81
## DF of independence model                     3.00
## RMSEA                                        0.00
## RMSEA lower 95% CI                           0.00
## RMSEA upper 95% CI                           0.00
## SRMR                                         0.00
## TLI                                          -Inf
## CFI                                          1.00
## AIC                                          0.00
## BIC                                          0.00
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)

## S matrix
mxEval(Smatrix, random2$mx.fit)

##           [,1]      [,2]      [,3]
## [1,] 0.7542121 0.0000000 0.0000000
## [2,] 0.0000000 1.0000000 0.1704927
## [3,] 0.0000000 0.1704927 1.0000000

## R2
mxEval(1-Smatrix, random2$mx.fit)[1,1]

## [1] 0.2457879

## Plot the model
plot(random2, color="green")

A path model for the cognitive ability to supervisor rating (1)

• The dataset is available in the metaSEM package as Hunter83.
• We may read the Excel data: Hunter83.xlsx.
• There are missing correlations in the dataset.

## Sheets including correlation matrices
sheet <- as.character(1:14)

my.R2 <- list()
for (i in sheet) {
  my.R2[[i]] <- readxl::read_excel("Hunter83.xlsx", sheet=i)
}

## Read study names and sample sizes in sheet "0"
my.study <- readxl::read_excel("Hunter83.xlsx", sheet="0")

my.R2 <- lapply(my.R2, function(x) {x <- unlist(x)
                                    x <- matrix(x, ncol=4)
                                    x <- vechs(x)
                                    vec2symMat(x, diag=FALSE)})

var.names <- c("Ability", "Knowledge", "Work sample", "Supervisor")
my.R2 <- lapply(my.R2, function(x) { dimnames(x) <- list(var.names, var.names); x})
names(my.R2) <- my.study$Study
my.R2[1:2]

## $`Campbell et al. (1973)`
##             Ability Knowledge Work sample Supervisor
## Ability        1.00      0.65        0.48       0.33
## Knowledge      0.65      1.00        0.52       0.40
## Work sample    0.48      0.52        1.00       0.23
## Supervisor     0.33      0.40        0.23       1.00
## 
## $`Corts et al. (1977)`
##             Ability Knowledge Work sample Supervisor
## Ability        1.00      0.53        0.50       0.03
## Knowledge      0.53      1.00        0.49       0.04
## Work sample    0.50      0.49        1.00       0.18
## Supervisor     0.03      0.04        0.18       1.00

A path model for cognitive ability to supervisor rating (2)

• Hunter (1983) conducted an MASEM to test a mediator model from Ability (A) to Supervisor rating (S) on 14 studies with a total sample size of 3975.
• There are two intermediate mediators, Job knowledge (J) and Work sample (W), in the model, whereas there is no direct effect from Ability to Supervisor rating.
• 

model <- "## Regression paths
          Job_knowledge ~ A2J*Ability
          Work_sample ~ A2W*Ability + J2W*Job_knowledge
          Supervisor ~ J2S*Job_knowledge + W2S*Work_sample

          ## Fix the variance of Ability at 1
          Ability ~~ 1*Ability

          ## Label the error variances of the dependent variables
          Job_knowledge ~~ Var_e_J*Job_knowledge
          Work_sample ~~ Var_e_W*Work_sample
          Supervisor ~~ Var_e_S*Supervisor"

plot(model, color="yellow", layout="spring")

RAM <- lavaan2RAM(model, obs.variables=c("Ability","Job_knowledge",
                  "Work_sample","Supervisor"))
RAM

## $A
##               Ability Job_knowledge Work_sample Supervisor
## Ability       "0"     "0"           "0"         "0"       
## Job_knowledge "0*A2J" "0"           "0"         "0"       
## Work_sample   "0*A2W" "0*J2W"       "0"         "0"       
## Supervisor    "0"     "0*J2S"       "0*W2S"     "0"       
## 
## $S
##               Ability Job_knowledge Work_sample Supervisor 
## Ability       "1"     "0"           "0"         "0"        
## Job_knowledge "0"     "0*Var_e_J"   "0"         "0"        
## Work_sample   "0"     "0"           "0*Var_e_W" "0"        
## Supervisor    "0"     "0"           "0"         "0*Var_e_S"
## 
## $F
##               Ability Job_knowledge Work_sample Supervisor
## Ability             1             0           0          0
## Job_knowledge       0             1           0          0
## Work_sample         0             0           1          0
## Supervisor          0             0           0          1
## 
## $M
##   Ability Job_knowledge Work_sample Supervisor
## 1       0             0           0          0

Random-effects TSSEM: Stage 1 Analysis

• There is a total of 6 correlation coefficients in the analysis with 21 elements in the variance component matrix of the random effects.

• We may specify RE.type="Diag" since there are only 14 studies.

• The minimum and the maximum \(I^2\) are 0.0000 and 0.7691, respectively.

• The random-effects model is more appropriate.

## method="REM": Random-effects model
random1 <- tssem1(Hunter83$data, Hunter83$n, method="REM", RE.type="Diag")
summary(random1)

## 
## Call:
## meta(y = ES, v = acovR, RE.constraints = Diag(paste0(RE.startvalues, 
##     "*Tau2_", 1:no.es, "_", 1:no.es)), RE.lbound = RE.lbound, 
##     I2 = I2, model.name = model.name, suppressWarnings = TRUE, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##               Estimate   Std.Error      lbound      ubound z value  Pr(>|z|)    
## Intercept1  5.0045e-01  2.8249e-02  4.4508e-01  5.5581e-01 17.7157 < 2.2e-16 ***
## Intercept2  4.3124e-01  2.1560e-02  3.8898e-01  4.7350e-01 20.0020 < 2.2e-16 ***
## Intercept3  2.0615e-01  2.5700e-02  1.5578e-01  2.5652e-01  8.0214 1.110e-15 ***
## Intercept4  5.2074e-01  3.1799e-02  4.5842e-01  5.8307e-01 16.3758 < 2.2e-16 ***
## Intercept5  2.5781e-01  3.5914e-02  1.8742e-01  3.2820e-01  7.1785 7.045e-13 ***
## Intercept6  2.3471e-01  1.7505e-02  2.0040e-01  2.6902e-01 13.4085 < 2.2e-16 ***
## Tau2_1_1    6.2569e-03  3.6590e-03 -9.1449e-04  1.3428e-02  1.7100   0.08726 .  
## Tau2_2_2    2.2234e-03  2.1500e-03 -1.9905e-03  6.4373e-03  1.0342   0.30106    
## Tau2_3_3    4.3637e-03  3.0130e-03 -1.5417e-03  1.0269e-02  1.4483   0.14754    
## Tau2_4_4    7.5512e-03  4.6477e-03 -1.5581e-03  1.6661e-02  1.6247   0.10422    
## Tau2_5_5    1.0692e-02  6.0658e-03 -1.1966e-03  2.2581e-02  1.7627   0.07795 .  
## Tau2_6_6    1.0005e-10  1.7130e-03 -3.3575e-03  3.3575e-03  0.0000   1.00000    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 205.225
## Degrees of freedom of the Q statistic: 60
## P value of the Q statistic: 0
## 
## Heterogeneity indices (based on the estimated Tau2):
##                              Estimate
## Intercept1: I2 (Q statistic)   0.7209
## Intercept2: I2 (Q statistic)   0.4475
## Intercept3: I2 (Q statistic)   0.5671
## Intercept4: I2 (Q statistic)   0.7459
## Intercept5: I2 (Q statistic)   0.7691
## Intercept6: I2 (Q statistic)   0.0000
## 
## Number of studies (or clusters): 14
## Number of observed statistics: 66
## Number of estimated parameters: 12
## Degrees of freedom: 54
## -2 log likelihood: -130.2043 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)

## Occasionally, we may encounter some error messages.
## We may try to remove the error messages by rerunning the model.
random1 <- rerun(random1)

## 
## Beginning initial fit attempt

## 
## Solution found

## 
## Solution found!  Final fit=-130.20428 (started at -130.20428)  (1 attempt(s): 1 valid, 0 errors)

summary(random1)

## 
## Call:
## meta(y = ES, v = acovR, RE.constraints = Diag(paste0(RE.startvalues, 
##     "*Tau2_", 1:no.es, "_", 1:no.es)), RE.lbound = RE.lbound, 
##     I2 = I2, model.name = model.name, suppressWarnings = TRUE, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##               Estimate   Std.Error      lbound      ubound z value  Pr(>|z|)    
## Intercept1  5.0045e-01  2.8249e-02  4.4508e-01  5.5581e-01 17.7157 < 2.2e-16 ***
## Intercept2  4.3124e-01  2.1560e-02  3.8898e-01  4.7350e-01 20.0020 < 2.2e-16 ***
## Intercept3  2.0615e-01  2.5700e-02  1.5578e-01  2.5652e-01  8.0214 1.110e-15 ***
## Intercept4  5.2074e-01  3.1799e-02  4.5842e-01  5.8307e-01 16.3758 < 2.2e-16 ***
## Intercept5  2.5781e-01  3.5914e-02  1.8742e-01  3.2820e-01  7.1785 7.045e-13 ***
## Intercept6  2.3471e-01  1.7505e-02  2.0040e-01  2.6902e-01 13.4085 < 2.2e-16 ***
## Tau2_1_1    6.2569e-03  3.6590e-03 -9.1449e-04  1.3428e-02  1.7100   0.08726 .  
## Tau2_2_2    2.2234e-03  2.1500e-03 -1.9905e-03  6.4373e-03  1.0342   0.30106    
## Tau2_3_3    4.3637e-03  3.0130e-03 -1.5417e-03  1.0269e-02  1.4483   0.14754    
## Tau2_4_4    7.5512e-03  4.6477e-03 -1.5581e-03  1.6661e-02  1.6247   0.10422    
## Tau2_5_5    1.0692e-02  6.0658e-03 -1.1966e-03  2.2581e-02  1.7627   0.07795 .  
## Tau2_6_6    1.0005e-10  1.7130e-03 -3.3575e-03  3.3575e-03  0.0000   1.00000    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 205.225
## Degrees of freedom of the Q statistic: 60
## P value of the Q statistic: 0
## 
## Heterogeneity indices (based on the estimated Tau2):
##                              Estimate
## Intercept1: I2 (Q statistic)   0.7209
## Intercept2: I2 (Q statistic)   0.4475
## Intercept3: I2 (Q statistic)   0.5671
## Intercept4: I2 (Q statistic)   0.7459
## Intercept5: I2 (Q statistic)   0.7691
## Intercept6: I2 (Q statistic)   0.0000
## 
## Number of studies (or clusters): 14
## Number of observed statistics: 66
## Number of estimated parameters: 12
## Degrees of freedom: 54
## -2 log likelihood: -130.2043 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)

## Extract the fixed-effects estimates
(est_fixed <- coef(random1, select="fixed"))

## Intercept1 Intercept2 Intercept3 Intercept4 Intercept5 Intercept6 
##  0.5004463  0.4312399  0.2061453  0.5207418  0.2578125  0.2347090

## Convert the estimated vector to a symmetrical matrix
## where the diagonals are fixed at 1 (for a correlation matrix)
vec2symMat(est_fixed, diag=FALSE)

##           [,1]      [,2]      [,3]      [,4]
## [1,] 1.0000000 0.5004463 0.4312399 0.2061453
## [2,] 0.5004463 1.0000000 0.5207418 0.2578125
## [3,] 0.4312399 0.5207418 1.0000000 0.2347090
## [4,] 0.2061453 0.2578125 0.2347090 1.0000000

Random-effects TSSEM: Stage 2 Analysis (1)

• When fitting a mediation model, we may want to calculate the indirect effect from Ability (A) to Supervisor rating (S):
  • Ind = A2J*J2S+A2J*J2W*W2S+A2W*W2S
  • We may request its 95% LBCI by specifying intervals.type="LB".
• The proposed model fits the data well with \(\chi^2(df=1, N=3975)=3.5748, p=0.0587\), CFI=\(0.9973\), RMSEA=\(0.0255\), and SRMR=\(0.0204\).
• The estimated parameters for the \(A\) and \(S\) matrices are and , respectively.
• The error variance in Supervisor rating is 0.9072. Thus, the \(R^2\) is 0.0928.
• The indirect effect (and its 95% LBCI) from Ability to Supervisor rating is \(0.1670\) \((0.1374, 0.1984)\).

random2 <- tssem2(random1, Amatrix=RAM$A, Smatrix=RAM$S, Fmatrix=RAM$F,
                  intervals.type="LB", 
                  mx.algebras=
                  list( ind=mxAlgebra(A2J*J2S+A2J*J2W*W2S+A2W*W2S, name="ind") )  )
summary(random2)

## 
## Call:
## wls(Cov = pooledS, aCov = aCov, n = tssem1.obj$total.n, Amatrix = Amatrix, 
##     Smatrix = Smatrix, Fmatrix = Fmatrix, diag.constraints = diag.constraints, 
##     cor.analysis = cor.analysis, intervals.type = intervals.type, 
##     mx.algebras = mx.algebras, model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: Likelihood-based statistic
## Coefficients:
##     Estimate Std.Error   lbound   ubound z value Pr(>|z|)
## A2J 0.510783        NA 0.456554 0.564987      NA       NA
## J2S 0.224033        NA 0.136802 0.311682      NA       NA
## W2S 0.121026        NA 0.056817 0.181657      NA       NA
## A2W 0.230913        NA 0.157756 0.299775      NA       NA
## J2W 0.397800        NA 0.311509 0.484687      NA       NA
## 
## mxAlgebras objects (and their 95% likelihood-based CIs):
##             lbound  Estimate    ubound
## ind[1,1] 0.1374292 0.1669702 0.1983626
## 
## Goodness-of-fit indices:
##                                                Value
## Sample size                                3975.0000
## Chi-square of target model                    3.5748
## DF of target model                            1.0000
## p value of target model                       0.0587
## Number of constraints imposed on "Smatrix"    0.0000
## DF manually adjusted                          0.0000
## Chi-square of independence model            975.2843
## DF of independence model                      6.0000
## RMSEA                                         0.0255
## RMSEA lower 95% CI                            0.0000
## RMSEA upper 95% CI                            0.0561
## SRMR                                          0.0204
## TLI                                           0.9841
## CFI                                           0.9973
## AIC                                           1.5748
## BIC                                          -4.7129
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)

## Display the A matrix
mxEval(Amatrix, random2$mx.fit)

##                 Ability Job_knowledge Work_sample Supervisor
## Ability       0.0000000     0.0000000   0.0000000          0
## Job_knowledge 0.5107833     0.0000000   0.0000000          0
## Work_sample   0.2309135     0.3978003   0.0000000          0
## Supervisor    0.0000000     0.2240334   0.1210258          0

## Display the S matrix
mxEval(Smatrix, random2$mx.fit)

##      [,1]      [,2]      [,3]     [,4]
## [1,]    1 0.0000000 0.0000000 0.000000
## [2,]    0 0.7391004 0.0000000 0.000000
## [3,]    0 0.0000000 0.6945954 0.000000
## [4,]    0 0.0000000 0.0000000 0.907194

Random-effects TSSEM: Stage 2 Analysis (2)

plot(random2, color="yellow", layout="spring")

Summary

• MASEM provides a useful framework to synthesizing and testing models in SEM.
• The TSSEM approach provides a statistically sound approach in doing MASEM.
• A new MASEM approach that handles continuous moderators is being developed. The paper and examples are available at PsyArXiv and Github.
• Another area of my interest is the SEM-based meta-analysis.
• Q&A!

---

1. National Research Council (1992). Combining information: Statistical issues and opportunities for research. Washington, D.C.: National Academy Press.↩

2. Brown, S. P., & Stayman, D. M. (1992). Antecedents and consequences of attitude toward the ad: A meta-analysis. Journal of Consumer Research, 19, 34-51.↩

3. Premack, S. L., & Hunter, J. E. (1988). Individual unionization decisions. Psychological Bulletin, 103, 223-234.↩

4. Norton, S., Cosco, T., Doyle, F., Done, J., & Sacker, A. (2013). The Hospital Anxiety and Depression Scale: A meta confirmatory factor analysis. Journal of Psychosomatic Research, 74(1), 74-81. http://doi.org/10.1016/j.jpsychores.2012.10.010↩

5. Murayama, K., & Elliot, A. J. (2012). The competition-performance relation: A meta-analytic review and test of the opposing processes model of competition and performance. Psychological Bulletin, 138(6), 1035-1070. http://doi.org/10.1037/a0028324↩

6. Viswesvaran, C., & Ones, D. S. (1995). Theory testing: Combining psychometric meta-analysis and structural equations modeling. Personnel Psychology, 48(4), 865-885.↩

7. Becker, B. J. (1992). Using results from replicated studies to estimate linear models. Journal of Educational Statistics, 17(4), 341-362.↩

8. Cheung, M. W.-L. (2014). Fixed- and random-effects meta-analytic structural equation modeling: Examples and analyses in R. Behavior Research Methods, 46(1), 29-40. http://doi.org/10.3758/s13428-013-0361-y↩

9. Cheung, M. W.-L., & Chan, W. (2005). Meta-analytic structural equation modeling: A two-stage approach. Psychological Methods, 10(1), 40-64. http://doi.org/10.1037/1082-989X.10.1.40↩

10. Cheung, M. W.-L. (2018). Some reflections on combining meta-analysis and structural equation modeling. Research Synthesis Methods, 0(0). https://doi.org/10.1002/jrsm.1321↩

11. Cudeck, R. (1989). Analysis of correlation matrices using covariance structure models. Psychological Bulletin, 105(2), 317-327.↩

12. Joreskog, K. G., & Sorbom, D. (1996). LISREL 8: A user-s reference guide. Chicago, IL: Scientific Software International, Inc.↩

13. Muthen, B., Kaplan, D., & Hollis, M. (1987). On structural equation modeling with data that are not missing completely at random. Psychometrika, 52(3), 431-462. http://doi.org/10.1007/BF02294365↩

14. Cheung, M. W.-L., & Chan, W. (2009). A two-stage approach to synthesizing covariance matrices in meta-analytic structural equation modeling. Structural Equation Modeling: A Multidisciplinary Journal, 16(1), 28-53. http://doi.org/10.1080/10705510802561295↩

15. Bentler, P. M., & Savalei, V. (2010). Analysis of correlation structures: Current status and open problems. In S. Kolenikov, D. Steinley, & L. Thombs (Eds.), Statistics in the Social Sciences (pp. 1-36). New Jersey: John Wiley & Sons, Inc.↩

16. Cheung, M. W.-L., & Chan, W. (2005). Classifying correlation matrices into relatively homogeneous subgroups: A cluster analytic approach. Educational and Psychological Measurement, 65(6), 954–979. https://doi.org/10.1177/0013164404273946↩

17. Jak, S., & Cheung, M. W.-L. (2018). Testing moderator hypotheses in meta-analytic structural equation modeling using subgroup analysis. Behavior Research Methods, 50(4), 1359–1373. https://doi.org/10.3758/s13428-018-1046-3↩

18. Rosenthal, R. (1991). Meta-analytic procedures for social research (Revised). Newbury Park, [Calif.]: Sage.↩

19. Schmidt, F. L., & Hunter, J. E. (2015). Methods of meta-analysis: Correcting error and bias in research findings (3rd ed.). Thousand Oaks, CA: Sage.↩

20. Michel, J. S., Viswesvaran, C., & Thomas, J. (2011). Conclusions from meta-analytic structural equation models generally do not change due to corrections for study artifacts. Research Synthesis Methods, 2(3), 174-187. http://doi.org/10.1002/jrsm.47↩

21. McArdle, J. J., & McDonald, R. P. (1984). Some algebraic properties of the Reticular Action Model for moment structures. British Journal of Mathematical and Statistical Psychology, 37(2), 234-251. http://doi.org/10.1111/j.2044-8317.1984.tb00802.x↩

22. Browne, M. W., & Cudeck, R. (1993). Alternative ways of assessing model fit. In K. A. Bollen & J. S. Long (Eds.), Testing Structural Equation Models (pp. 136-162). Newbury Park, CA: Sage.↩

23. Digman, J. M. (1997). Higher-order factors of the Big Five. Journal of Personality and Social Psychology, 73(6), 1246-1256. http://doi.org/10.1037/0022-3514.73.6.1246↩

24. Epskamp, S (2017). semPlot: Path diagrams and visual analysis of various SEM packages’ output. R package version 1.1.0. http://CRAN.R-project.org/package=semPlot↩

25. Becker, B. J., & Schram, C. M. (1994). Examining explanatory models through research synthesis. In H. Cooper & L. V. Hedges (Eds.), The handbook of research synthesis (pp. 357-381). New York: Russell Sage Foundation.↩